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ABSTRACT 


A new method of numerical solution for time-harmonic 
field scattering problems involving inhomogeneous bodies of 
revolution is presented. A simple single-loop multiport 
feedback system is driven by the incident field. The forward 
path transfer matrix, found by use of the finite-element 
method, translates the total near-field into the body-surface 
field. The feedback path matrix is established by the magne- 
tic vector potential formulation; it translates the body- 
surface currents into the scattered near-field. The thin- 
wire scattering problem is formulated in terms of this field 
feedback method. Results of this computation are compared in 
graphical form to those found by Hallen's integral equation 
solution of the scattering problem for several cases of 


ipeerest . 


E SCHOOL 








i 


ite Or CONTENTS 


IT. INTRODUCTION- - - - - - ------ - 


eer ROBLEM AND METHOD OF SQOLUTION- - - =---- - 


eer ini D FTEepBACK FORMULATION- - - = -- - - 


Pepe, WIRE, SCATTERING PROBLEM - 


Ge. COUPLED AZIMUTHAL POTENTIALS- - - 


me mm mm mm i i SE Sl i= 


III. COMPUTATIONAL PROCEDURES- - - - - -------- 


ee Peet INTE SEU EMENT METHOD - 


1. Basis Function Expansions - - 
2. The Euler-Lagrange Variational 
Technique - -----+e---- 


3. Application of the Thin-Wire Problem- - 


Pees SCHIiEReED FIELD INTEGRAL 


L. The Vector Potential Formulation- - - - - 


2. The Scattered Field Matrix Equation - - - 


C. PROBLEM IMPLEMENTATION 


TV. RESULTS AND ANALYSIS- - - - - - - - - ---- we 


A. INTEGRAL EQUATION COMPARISON 


B. ERROR ANALYSIS- - - - - - = ---+-e-e-e-e- - 


V. CONCLUSIONS - - - --- - 


APPENDIX A- - - - =---+e-ee-ee-ee-8e8e8-8-28-4e82-8282-8- - 


APPENDIX B- - - - - ----+-+--e-4--e8--4-4-8- = 


APPENDIX C- - - - - - - - = = = 


APPENDIX D- - - =- = -=--+-+-e-4--e-8--8-e8-4-428-8- 6 


LIST OF REFERENCES- - - - = - = = = = -------2e- - 


INITIAL DISTRIBUTION LIST - - - 


10 


eZ 


ne 


16 


i 3 


oe 


24 


ae 


oy 


Ze 


a3 


38 


39 


U3 


47 


ug 


87 


105 


odes 


1338 


136 





ioe ln PRODUCTION 


The objective of this endeavor is to investigate the per- 
Formance of a new, and potentially quite powerful, technique 
fOr the computation of electromagnetic scattering problems. 

For reasons that will later become apparent, this new technique 
feeeebe denoted the "rield feedback formulation" (FFF or on 
The work reported here is a first application of en where 

the calculation of scattering from straight thin-wires has 

been performed. This is an almost classical electromagnetic 
field problem that has been routinely solved by use of inte- 
gral equations for many years. Such a cononical problem was 
chosen as a testing vehicle for the newly conceived es 
analysis method since alternate integral equation solutions 
are readily available for comparisons. 

The formulation developed here employs a combination of 
the finite element numerical method and the classical magnetic 
vector potential. The finite element method (FEM) 1s an 
application of variational techniques which implements a 
weighted residual approximation. The high accuracy of the 
finite element method is used to advantage in reducing the 
amount of computing required, as compared to other numerical 
methods. Use is also made of coupled azimuthal potentials 


(CAP). The CAP formulation restricts the class of approachable 





problems to those which are axially symmetric, but does allow 
the treatment of material inhomogeneities and multiwavelength 
dimensions. 

Thin-wire assumptions are made which greatly simplify 
the current-scattering problem. These assumptions, however, 
are by no means critical to the solution method. The formu- 
lation should be readily extendable to a large class of impor- 


tant interior, radiation, and scattering problems. 





eee ieee eoouenh AND METHOD OF SOLUTION 


Flectromagnetic scattering problems usually involve 
a known incident field, such as a uniform plane wave. The 
scattered field can be found from the current induced by the 
presence of the incident field, but the induced current can 
be found only after the total field is known. Since the 
total field is just the sum of the incident and scattered 
fields, a definite difficulty has been encountered. Know- 
ledge of the induced current is needed to find the scattered 
field, while at the same time, the scattered field must be 
known to determine the induced current. 

This is clearly seen in the form of Maxwell's equations 
for the scattering problem. Assume that a scatterer with 


properties Y +3 WE Z,=jwu, is in free space where 


im 1 lee ue 


ee OU. > ZF JOU,» and in the presence of sources Jo» M 


=O 


O° 


Then the curl equations have the form 


Peo =o) + Y_ © (la) 
e) IE 
VxE = -M -A,H Gia) 
or il. 
where the induced electric and magnetic currents are 
Ji = (Y,-Y QE (2a) 
My = (2, -2.)H 





and the scattered magnetic field is 


— —x; — 


= H- 4H Cay) 


H , J 
scat incident 


Miewelectric field relationship is, of course, similar. 

The simple and compact nature of the simultaneous equa- 
femens (1) and (2) is highly appealing, but of little practi- 
cal use in obtaining numerical answers to a problem because 
of the excessive computing required. Yet, while the require- 
ment to have two unknowns cannot be escaped. the equations 
relating these unknowns can at least be rendered more 


tractable. 


Pees FLELD FEEDBACK FORMULATION 

As will weston later, two simple matrix equations can 
be derived which relate the magnetic and electric fields at 
a number of points along the surface of the scatterer to the 
fields at points along a boundary some distance away. These 
equations are independent and work in opposite directions. 
Tmet as, from the first equation, given the total boundary 


Gmetad values, e, and Ay» the scatterer surface total field 


c 


values, e andh_, can be found: 
—W —) 
; ee | a 
Lew h ] - V fe, hd (4a) 


I 


The column vectors, om. contain the field values at the 
SiS 


computation points. The second equation takes fe ih, JP and 


yields the scattered field values fe_'h]” along the boundary: 
| 





J (4b) 


Equations (4a), (4b), and (3) comprise a simple feedback 


Pyetem, as in Figure 1. 





Figure 1. Field Feedback 


What remains, then,is to construct the matrices V and T 
memmete tne problem at hand. It should be pointed out that 
the magnetic field alone 15S involved in the wire-scatterer 
problem as developed here. So Ceth J" is replaced by hy: 
The field translation matrices V and T, though developed in 


exactly the same way, will be considerably simplified as 


Compared to the general scatterer case. 


Peelot WIRE SCATTERING PROBLEM 
Given a perfectly conducting thin wire in free space, and 
a plane wave of known amplitude incident from a known direction, 


it 1s required to find the current induced on the wire and the 


10 





Pe wamaecaeteredq by the wire at any point in space. Finding 

mpemenaucead Current 1s equivalent to finding the field on the 
A A 

maeeensiuntace, since J = nx H , where n is the unit vector 


Memmel tO the surface. Normalized coordinates will be used 


throughout the problem, and are given by 
mon.) = (Kor, K G2 o>). 


The coordinates (r,Z,¢) are standard circular cylindri- 
cal coordinates, and K Fema. is the free space wavenumber of 
the time-harmonic field. The situation is depicted in 
Figure 2. Use of these normalized coordinates will allow 
the formation of a compact and useful form of the CAP equations. 

The wire radius 1S a, given in terms of ae the free space 
wavelength. Then the normalized radius is Ri=27a/h_. The 
wire length is L, also given in wavelengths. So Z,s2mL/r_. 
The incidence angle of the plane wave is a, measured as 
shown in Figure 2. Since perfect conductivity is assumed, 
meas and Currents will exist only on the surface of the 
(eer end Since a very thin wire will be assumed, where a 
<A /10, the current at the wire ends may be set to zero. 

The problem will not be solved in terms of the usual 
electric and magnetic fields, E and H, or the scalar poten- 
tial, V. Instead, use will be made of the magnetic vector 
potential, A. The current and scattered fields will be 
presented and calculated in terms of the CAP formulation 


discussed below. 


ile 








Figure 2. Normalized Coordinates 


SeeecOoUPLED AZIMUTHAL POTENTIALS 

The CAP formulation [3] has been developed to handle 
previously intractable time-harmonic field problems involving 
continuously and discretely inhomogeneous bodies of revolu- 
tion. With this limitation of axial symmetry, all electro- 
magnetic field components can be expanded via exponential 
Fourier series in the ? coordinate. The model fields can 
then be represented by two coupled azimuthal potentials. 
These potentials satisfy a system of coupled partial differ- 
ential equations, as well as a variational criterion, which 
Will be developed here. 

Using the normalized coordinates in Figure 2, the total 


fields are decomposed into azimuthal modes: 


ey 





more Zo) = enn (R,Z)exp(jmo) CS) 


m=-—% 


Nn H(R,Z,9) Z n,(R»2exp(jmd) (5b) 
where nywi2er. By substitution of these expansions into 
Maxwell's equations, it can be shown that all modal field 
components can be generated from two coupled azimuthal 


potentials ), (R,Z,m) and Vo(R,Z,m), using 


4 x en = j£,,(moxVb,-Ru Vy) (6a) 
$° Es v,/R (Sb) 
4 x he = j£,,(moxVp,+Re_ Vp.) (Ge) 
, h = Wo/R (6d) 


where ¢ is the unit azimuthal vector, V is the two dimen- 


Slonal gradient operator, and 


: Zea — i 
FUR 2Jd=lu (R,Z)e (CR, Z)R -m° J (7) 


The Euler-Lagrange variational technique will be dis- 
cussed more fully in the exposition of the finite element 
method below. It is sufficient to note here that this tech- 


nique finds the solution of the system partial differential 


des 





Pa@eetona as the Stationary point of a functional of the 
potentials and their first derivatives. This stationary 


point is defined as the point (Ws 9: ) at which the 


Yoo 
variation of F, 6F, vanishes for perturbations (Op, Sw.) 
Beouwesthat point. 

The functional F is a surface integral of the Lagrangian 


L, where 


Ee = £ WLVp, (Re Vp, tmoxVy, +0, * (RUVo,-moxVy, DJ 


Zz Zz 
- Ce vy tL, YR (3) 


The functional F is integrated over the constant azimuth 
planar meridian cross section 2 of the volume in which the solu- 
imemeis Of interest. 

It 1S convenient at this point to change variables, and 
introduce a simplifying assumption. From (6b) and (6d) we 


have e a R222 =H4/R and h at R24) - Wo/R. So) le sjoleliegless epnliy’ Telave 


p p 
azimuthal components of the modal fields, atm and Dom? will 
yield a complete solution through (6a) and (6c). Now, since 


a thin wire is assumed, the phase variation of the incident 
field from side to side on the wire may be neglected. Then 
miemimncuced current will be a function of Z alone; cee 
Therefore the scattered fields will be azimuthally invariant, 
pememat ECR,Z,0) = e(R,2) and nN H(R,2,¢)=n (R,2). Further- 


~ 
heme, samee J_=nxH=H,Z and since 


Z Q 


14 





= alk = ilk 2 dp a el ee: 
JWE oe JWE L=R 0Z R oR R Hyd, 
fam Olilows that e Swiemeanad h = h . The number of modes 
bo O oOo 
has now been reduced to one, the h mode. 


dO 


Observe that (6a) and (6c) no longer involve coupled 
potentials. This allows the Lagrangian to be decoupled, and 


it may now be written 


a 2 a 2 
—- £,tV(Rh, ,) ERY CRAY) It MARR Go: (9) 
Bowene functional is 
= f LCR,Z,h, Why dRdzZ 0; ) 


22 


memeemains to find Ngo Rs2) EOunicncCe wesc aALIOnNary in {) for 
given Dirichlet boundary conditions. This will be done by 


use of the finite element method. 


15 





TIt. COMPUTATIONAL PROCEDURES 


mewcoustructions of the VY and IT matrices have little 
in common other than the coordinate system. The finite 
element method is used to generate V. T 1s derived from 
the magnetic vector potential in this perfect conductor case, 
while if a dielectric body were used, both electric and 
magnetic vector potentials would be required to construct tT. 
The use of numerical methods obviously implies that the 
current and fields will be found only at discrete points. 
Pemeceseussed in Section I1.C., only points in a constant 
azimuth meridian plane need be considered. The region between 
the wire surface and the boundary at which the fields are to 
be calculated is divided into triangular elements. These are 
areas of integration in the finite element method. The divi- 
Sions and the nodal numbering scheme are shown in Figure 3. 


There are a total of Nv Pewee, alone the wire, and N. points 


B 
along the boundary. 


Pee UNITE ELEMENT METHOD 

This method is one whereby the region of interest is divided 
into a number of overlapping subregions of a specified shape. 
Then the problem solution function § is approximated by 
defining the solution as a linear combination of the products 


MEmiienown coerficients with "basis" funetions defined over 


16 





—————— 


over each of the overlapping subregions. These basis 
functions are usually sectionally linear, and, in very simple 
terms, are used to "fill in" the approximation between points. 
A weighted residual approximation of the solution 9, commonly 
known as the Ritz method, can then be made by the application 
of the discretized Euler-Lagrange variational procedure, 
subject to the appropriate boundary conditions. The rigorous 
development of each of these parts of the finite element 
method is not necessary here, but the background required for 


a basic understanding of the method will be given below. 





Pastimes sae SOlutoneroimtrs 


ley 





esas hinecrion Expansions 

imperder CO solve a problem numerically, some method 
of approximating the solution using a finite number of 
Mmenowns must be found. Also, the method's convenience must 
be baianced against the desired accuracy. A general term 
aeemsweh an approximation is a "basis function expansion" 
Semewitunection 6 |2]. The type of expansion used in the 
finite element method for a one-dimensional problem is given 


by 
N 


8(X) = > 8a, (X) (11) 


1=1 


where oy is the discrete version of 9, the os are the 

values of 6@(X) at the discrete points Xe and the U,¢X) 

are the basis functions at the points x. Linear basis 
functions are usually employed in FEM calculations, though 
quadratic or higher-order polynomial functions could be 
used to obtain greater accuracy. Linear basis functions for 
a one-dimensional problem are shown in Figure 4. Note that 
the node spacing is irregular. This allows tailoring of the 
expansion to better approximate the solution, and, in two 

or higher dimensional cases, may be necessary because of the 
shape of the solution region. If neither of these are 
necessary, regular node spacing can be used with some advan- 


tage to reduce complexity. The basis fucntions u, are seen 


18 





to have unit value at their respective nodes Xa and are 


zero for xX. 


yo xX. 
1-l- — 


esteem Li <N 


a 


Figure 4. Linear Basis Functions 


A two-dimensional linear basis function often used is 
the "pyramid" function; its two-dimensional projection is 
shown in Figure 5. A rectangular region 2 in the x-y plane 
has been discretized by division into the triangular elements 
shown. The boundary 82.3 5 wise AV) oumemeclal=: Oaibhalens slong U3 96%) has 
been outlined to make the shape clear. 

The function has a value of one at the center node, and 
is zero along and outside the periphery. Its equation can 


be written 


Uz, (x,y) = oe Cre) Meals Ci? ) 
where the component functions are planar segments of the form 


ue Meme tS 7 ty Cir) 
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Rectangular 
region 2 





3 
2 
ie L 
7 ey. 3 uy 5 6 7 x 
1=1 
Figure 5. Two-Dimensional Finite Element Mesh 


Peeeces ,y must be functions of (x,y) so that u2 


396%30V_)1, 


ud 9 (xs 9y5)=0, ef o4 it?) Bem intecen (i 4) 

As in the one-dimensional case, the basis functions over- 
lap; that is, a function is centered around each node in the 
region or on the boundary. Of course, on the boundary of 2 
meme dae half of the function is "lost" outside 2. The 
overall appearance of the basis function array in three di- 
mensions might be imagined as a crystalline surface of multi- 
faceted peaks centered over the nodes, with raised valleys 


running between. If a function @(x,y) were defined over Q, 


20 





its expansion could be written 


a 6 
9 6%2¥) = 2, ae en ce (x,y) (14) 


i oe 
If the exact values of 954 were known, 4 would be a multi- 
faceted approximation to 3 » egdualk only at the nodes x... 
a) ae 


Then, within the limitations of machine accuracy and computing 
time, the greater the node density the better the approxima- 
een O45 ei be. 

The wire scattering problem expands the field, Noe? by 
use of these pyramid basis functions. Referring to Figure 3, 
it is seen that the region 2 covered by the triangular ele- 
ments will not contain any complete pyramids, but this is not 
meeem@eyemecessary. it will prove convenient to use a sequen- 
tial node coordinate system rather than the (1,3) Cartesian 
description of (12). Let n be the node coordinate, where 
l<n<N, N= =NtNa- Then the approximating expansion Don of NG 
can be written 


N 
G4) = of eee CR AZ) Cis) 
non 
n=1 


Don 


Note that the top and bottom wire nodes are not represented. 


This is because of the assumption that h,_=0 at these points. 


dO 
Two field vectors can now be defined. The nodal wire- 
surface field is 
h =h "ane = i. N (16) 
—w —n W 


Za 





The nodal boundary field is 
ae = nh » ne H tl, N Cs) 


These vectors contain the approximate values of h at the 


oO 


nodes of the triangular mesh. 
feeesuler-Lagreange Variational Technique 
Consider the differential equation A9=f, which together 
Paeeeertain boundary conditions, describes some physical pro- 
blem in a region 2. The complex scalar wave (Helmholz) 
equation, A =[y +k’ ]6=0 1S an example. The differential equa- 


tion is related to the functional F(8), where 


F(@) = J “(L(8, DS, D8, x, y)dn Gis) 
42 

This relationship is a fundamental result of the calculus of 
variations, and simply stated, is the fact that the functional 
Mapes Stationary at the "point" 8 which also satisfies the 
Seeeperemciad! equation. That is, if the particular function 8 A 
can be found which causes the derivative (or first variation) 
ees to vanish, then FC@) is sadd to be stationary, and 
8S mcd to be the stationary point of F(6); ae is then 
ie Selution of the equation Aé=f {({5]. 

Observe that the vanishing of dF/d@ implies that a 
Pana, Maximum, or saddle point of F occurs ae In many 
physical problems F is associated with the potential energy 


of systems and thus finding oe 1s sometimes termed the 


Z2 





"Minimization" of F. The problem now presents a choice of 
computational methods. The differential equation might be 
approximated by a discrete system, say finite differences, 
er the variational integral, F, can be "minimized" over a 
finite number of trial functions, as is done by the finite 
element method. 

It is necessary in this procedure to find the Lagrangian 
Peeeoreiise in F(6) [4]. This L must be such that its insertion 
ime the eee eeGnes Sauaition produces the differential 
equation that is to be solved. As an example, the Euler- 
Lagrange equation for the two independent variable case is 
given by 


OL 3 rok 


59 7 ON D8" 0 (19) 


eg co 

oL D,,§ 
There is no definite way to discover L for any given problem, 
but it usually yields to informed guessing. 

The problem's boundary conditions present no unmanageable 
Ceeereutle2es in this formulation. If 6 is specified on the 
boundary (Dirichlet boundary conditions) nothing else need 
be done. If 9 is unspecified on the boundard, the stationary 
solution will automatically fulfill the Neumann boundary 
eomda tion, go =0, also known as the natural boundary condition. 

dn 


Boundary conditions of a more complicated type will require 


the functional F to be modified, but this procedure will not 
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Pee eescery mere, rE can be easily shown by using the field 
generating equations in [6], that the natural boundary condi- 
tion here results in a zero tangential electric field on the 
wire surface. 
3. Application to the Thin-Wire Problem 

pemdiscussed in Section II, the driving force in the 
scattering problem is the incident plane wave. Then, the 
boundary conditions to be applied are the natural condition 
along the wire, for nodes n=l, Noy dviceed Dirichlet "ond meron 


along the outer boundary of 2; h atte Creal field). Uae 


bo 
Peon lLet Condition Dee has also been imposed on the top 
and bottom wire nodes. These correspond to boundary condi- 


tions along the top and bottom node rows of a larger mesh, 
as in Figure 5. 

The minimization of F for the wire scattering problem 
mem@adrried Out as follows. F is discretized by substituting 


Don fOr Deo? 


which yields the approximate functional 


c c . 2 
FChgy d= f CER? RU V CRAG) F-Ru hy IdRdZ @20 ) 
2 


iim= exact functional F(h Be is made stationary about the 


9 


solution Ne by differentiating F with respect to hy, Se las 


makes F stationary about every point on Neo Since the appro- 
fumes Solution is to be found only at a finite number (N_) 

wl 
of points, the functional needs to be made stationary only at 


those points. This is done by a simultaneous solution of the 


No stationary condition equations 
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dRdZ, i =N om 
W 


where he 1s the unknown value of h at the node n=1, and 


ON 


Wal. Writing h explicitly, 


oN 
N N 
a . } ; ) - a2 
fem; VR hu -Ru,; hu, lJdRdZ=0, 1=1,Nw (22) 
2 n=1 n=l 
N 
Gua) since ) he is a linear operator with no coordinate 
n= 
Sep-nacnce, 
N 
> h y [=VRu.-VRu_-Ru.u_]dRdZ=0, i=1,Nw Ge) 
a R a n oe 
n=1 02 


Boundary conditions at the top and bottom wire nodes 
are already included in these equations. Application of the 
meet telad boundary condition is accomplished by observing 
that the ho» atl SN are known, where, using (15), hy=h,. 
So by splitting the sum the useful minimization equations 


for F (nh N? are found to be 


® 
N N 
Val 
@, he > C.h = 0 , i=l,Nw (2ha) 
ie D Be at 
= n=N +1 
W 
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where 


Cs iy, [=vRu.-VRu -Ru.u_JdRdzZ (24b) 
iT) R is n ies 
Q 


The matrix equation to be used for computation can be 
memmed directly from (24). First, recall that the basis 
functions u, are zero on and outside of the boundary 322 


This means that the product of non-overlapping basis func- 


mieeemis Zero. For example, in Figure 3, uz°u,=0. Thus 
mee = 0 , jjn-i| > 2 (25a) 
1, Vie! — 
and 
VRu.:VRu_=0 , |n-i| > 2 (28) 
sf n ~ 


This behavior is very useful in that it produces banded 
coefficient matrices. Special techniques can be used to 
invert large banded matrices very quickly. Writing out a 


few equations should make the matrix construction clear. 


fee eco nec on +c "7 : Oe (0 
ie) ie IN +1 ON +1 IN +2°N +2 
meade 1 2097 “pata * “on +2°y +27Son tay +3 = 9 
W W W Ww 
= > ; 
L=32 0 CoohytCg anata hy + Con a3hy 43°C3n 4uPy +4 = 9 
W W W W 
(26) 
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The system equation can therefore be written 


ime 2D, a? 
where the NAN GGcritelemm matrix A contains the basis 
momection integrals for 1=1,N ,n=1,N and the N XN TL matrix 
fee om@tains the integrals for 1=1,N,nsN +1,N. Both A and B 


are diagonal and have the forms 
le ine 


X Xx XX 
A: ioe 
y a 0 x a 0 
x xX X cae 
me, xX Xx | er 
x x ane 


The wire field ie due to the boundary field h, 1s now 


known, since 


ae ic : 
h =-A “Bh, = V h, (28) 
Equation (28) is the forward path in the system of Figure l. 


The actual calculation of the Asn Bs Will be discussed in 


Seemionm LLI.C. 


eeoece TERED FIELD INTEGRAL 
1. Vector Potential Formulation 
The development of the feedback path is less compli- 


@aeee-c eam the forward path construction given in III.A. The 


PE 





thin-wire approximation is applied to the familiar classical 
vector potential A, and a matrix equation formulation built 
on a one-dimensional basis function expansion. The result 
will be an approximation giving the near field produced by a 
wire-surface field or current. 


From the inhomogeneous Helmholz equation 
VAt+ kA = -J (28) 


mes seen, since JeJ5Z by assumption, that A=A_Z. Then it 


Z 
follows from the vector potential defining equation 


fee Vv x A Go) 
that 

pees A (R.Z) (31) 

““¢ 7 OR a 


The azimuthal invariance of H has been made explicit in (31). 

A commonly used expression for A [1] is given as an 
integration of a weighted point-source Green's function G(r',r) 
Over the current source where 


eli. | aan 


-, —-, _ ll 
aoe ae) = at 7 eieerl: G32) 


Ga) = | Tine Giant a) dy" (33) 
t 


v7 


This general situation is shown in Figure 6. 
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Figure 6 


2. Scattered Field Matrix Equation 
The formulation is easily applied to the thin-wire 


case. Recalling that h = no , and assuming that the 


pO ® 


eweeent J5 lies along the Z axis, Sno unom thin wire. ss 
; Di, was 12 | 
Nofz (RZ) = hy &Z ) ort R_dodd (34) 
S 


where, ID] = vR7+(Z-2) 4 ; as seen in Figure 7, and the integra- 
Yen 1S over S, the wire surface. Now, putting (34) in (31) 
yields 


RR fl/2Z 
eGR «2 2——— ae 
O 2 


p 
-1/2 ZF 





Ca J ia 1 4673 1Pl az: C35) 


é eee) 
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Figure 7 


This Ryo SR>2) 1S approximately the scattered 
eae , nh(R,2), radiated by the wire carrying an induced 
current Jo: Wits Wire Clewent, in Ene form of eee must 
be approximated by a basis function expansion so that nh (R,2Z) 
can be computed. Since ey oe Tose eUnetion) OL. 24 On lynne 
one-dimensional linear basis functions shown in Figure 4 may 
be used. Then, for the same ie nodes along the wire, as in 


Beerion LIt.A., 


5 a") = ez) (36) 
O EES 


p 


fimemee that (36) defines the same nodal surface-field vector 


eweelo). Substitution of (36) into (35) produces 


oo 


| 








! 
RR Ma antl 
pmGR. Zz) =—-— h eC ees je IP lag: 
S 2 n a Zz eS 
= 
n Ze mes 


|D| | D| 
Ga) 
ues . : 
ime limits of integration reduce from +2, 79 to 2 ae, 
1 ; 7 : ! . cs ? ' 
Z etl because the basis function u_oZ eels eZero “few Zt 124 2Z a1: 


Equation (37) can be used to construct the T matrix by noting 
tmat (3/7) gives the scattered field at one point (R,Z). Then, 
by choosing a number of points in some region, a discrete 
field approximation can be constructed. The points chosen 


Will be, unsurprisingly, just those Ne points along the finite 








element mesh border. So the matrix equation is 
eee i (38a) 
—S ~ —W 
where 
Zo 
een : | 
c= =| UH eae =] eI IPl az: (38b) 
' |D| |) 
ni=- 


is the coefficient of the n-th value of the wire field, and 


ip Memene term Of the b=-th value of the scattered field. 


Bron 


meneal ily, Since nh = h_ + Maes 


o>? an expression for the 


Eesti ela in terms of the incident and wire fields is 


h, = h, ae ee) 
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The scattered field ie found in (38) is a near field 
expression. It is valid at any distance from the wire, but 
this behavior does not extend to the entire problem. If he 
were to be found at some large distance from the wire, simply 
extending the radial width of 2 would not suffice. This is 
because FEM calculation accuracy degrades as the distance 


from the wire to the boundary, B increases. The magnitudes 


1? 
@r the electrostatic and induction fields decrease rapidly 
as R increases, so the point is quickly reached where they 
Beeome effectively zero as far as computer calculations are 
concerned. The resulting incorrect total field values are a 
part of the FEM R-related inaccuracies. The greater part of 
GeenrEM error is found in the basis function expansion of Ny ° 
The planar-section basis functions used are unable to 
accurately model the near field over a large radial distance, 
where "large", in this case, may be 1/100 of the wavelength. 

To solve the far-field problem, the near-field solu- 
tion must first be obtained for the wire field Le Themes 3) 
may be used to find he along any boundary. For typical far- 
field problems, assumptions may be made to reduce the integral's 
complexity. Let 8B be the angle of the far-field point, 
measured from the Z axis at the origin. Then R=kor Sire , 
where r is the spherical radius of the point. The far-field 
assumption is 


|B kor - 2! cosé (in phase term of (38b)) (40a) 


|D | 


KGP (in denominators) Cup) 
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50 the radiation field along a spherical are of radius r is 


given by (38), where Thin 1s 
ee 
mer 
Med Still. aicer Peni Oconee 
e ee O | a2 de dZ Ga) 
t 
7 n-Ll 


Seen ObLEM IMPLEMENTATION 


Combining equations (27) and (39) gives 


ia = ~(A+BT) 1B 


ie, Bl C42) 


on 


ieee 


This is, of course, the same equation as is represented in 
-1 


Figure 1, where V = -A “B. The system "transfer matrix" 
is defined by 
meee y Ty Vv eh (43a) 
ew Ee 3 % ~ = rie 
Simple matrix algebra leads to a simplification of this 
eesult 
me BT) CB Ch. (43D) 
—W ~ ~ ~ ~ ale ae 


Multiplying this equation by (A+BT) yields (42), so the system 
equation as derived does behave as a feedback system. 


, the 


The only variable remaining to be identified is ae 


incident field vector. The incident plane wave Bis can be 
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expressed in the normalized cylindrical coordinate system as 


h,(R,Z) = je 140080 5 1 (Reina) (4) 


where Je! is the derivative of the zeroth-order Bessel func- 
[em olnce R will be kept small, as discussed in Section 
TIT.B., (44) may be simplified by replacing Te! (Rsina) 


with the first few terms of its series expansion 








e e e 5 
: : Rsina Re ehO see stan 
t a 
Then Ln Mayo be cCalculared at the boundary nodes by use of 


C44) and (45). 
Now that pe eset One;EOD ii vOGwchne Guantiaities in G42) have 
Peemeround, the computer program can be devised to carry out 
the calcuations. The use of an existing routine, CSMINV, for 
the inversion of (A+BT) reduces the number of major computations 
to two: equation (24b) for the A and B elements, and equation 
(38b) for the scattered field coefficients Ton: 
The integration of (34b) becomes alzgebraicly messy, but is 
otherwide straightforward. Because of the basis function's 
behavior, the integration area is reduced from 2 to 2, the 
region of non-zero overlap of us and UL. This integration 
over o2 may then be written as a sum of integrations over the 


triangular elements comprising ®, since ue ls plecewise planar 


in the elements. Also, note that or LS independent of Z in 
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SS eS esr en eee 


the sense that the integration for any pair (i,n) of nodes, 
will be the same if that pair is translated along the wire. 
This means that one wire node may be selected, and the inte- 
grails concerning that node and its immediate neighbors will 
be the only ones required to construct A and B. For example, 


Beeect 1=2. Then C € c and C are the only 


ieee 2 2 ZN ie 
W W 


integrals that need to be calculated. (Node three is not 


being neglected. Equation (24b) shows that C = C leibie 


Zo ae 


mime its just C moved up one element. SoC C ) 


21 ome 


In this simple narrow mesh case, the integration of (24b) 


1s required over at most three elements (C,,), and the size 


22 7 
and shape of all elements is the same. For many problems 
this will not be the case. For large meshes, and especially 


for skewed elements, it becomes very convenient computa- 
tionally to change from a global to a local node system. 
Instead of moving from node to node, calculating all the 
coefficients about each node one at a time, an algorithm 
is devised to move from element to element, finding the 
integrals involving the three corner nodes over each element. 
This procedure effects a considerable simplification in the 
“bookkeeping” necessary to keep track of the nodal indices 
and coordinates. 

The local node system employed here is shown in Figure 8. 
The integrals of the combinations of local nodes (k,l) are 
found over a upper and lower element. Then the integrals 


can be placed in their proper global positions. As an example, 


ao 





Upper 


Element 
Lower 
Element 
Figure 8. Local Node Indices 


consider the upper element with global nodes C25 ee The 
integrations of (24b) over U are found for local node pairs 
Memes 2),(1,3),¢03,1),¢03,2), and (3,3). Then the local 


Me@embpair (1,1) integral, f(1,1) is put in A as one part in 


three of Cy 9: POU) ee smote Dart 1 Eworor Con, +3 in B, and 
f(1,3)=C,,. aT A. Also, f(3,1)=C,,, Gna) ols part OF CN, +3? 


meee (3,3) 1s part of C The remaining three local node 


oon 
integrals over U are not needed in this case, but would be 
used in larger meshes. 

For the thin mesh used here, the local node procedure 
improves only slightly the algorithm efficiency as compared 
to a global node procedure. The principle reason for its use 
was that an existing and successful subroutine was available 
to do local node integrations. This routine is named VARELA 


immemne program. Jt requires only the local node ordering 


scheme and the coordinate (R,Z) of each local node to produce 
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a three by three matrix containing the integrations. As 
discussed above, the wire-scattering problem needs only two 
calls to VARELA to obtain the necessary values. The sub- 
routine ABLOAD then fills A and B with these integrals. 

The calculation of the coefficients Ds is done by the 
familiar Gaussian quadrature method, using eight unequally 
spaced points. The interval of integration in (28b) is split 
into two segments, as in Figure 9. Then a change of variable 


is affected so that the limits of integration in each sub- 


interval are X=+l to -l. The Gaussian procedure then gives 
pe, n 

| f(x)dx = wiLi(x,)+f(-x,) | (46) 
—1 ee 1 


where the values tx. and the weight factors Ww. are tabulated 

in several references such as [6]. The code used to imple- 

ment these integrations is contained in subroutine GTrINT. 
moeedon (42) is coded in two subroutines, CFORM end 

ZCURR. Routine CFORM produces the matrix product 

(A+BT) 7 "B=C, while ZCURR calculates the incident field 

vector ae according to (44) and (45), then finds the wire 


fmrera vyector h =C h. : 
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IV. RESULTS AND ANALYSIS 


In many other time-harmonic field problems, FEM calcula- 
tions yleld good results with perhaps 20 to 30 nodes per 
wavelength. With higher node densities, say 50-75 perwavelength, 
the solution breaks down due to roundoff error in the compu- 
tations. As will be shown, this does not appear to be the 
case with the field feedbaek formulation. Nodal densities 
as high as 100 per wavelength have been used without diffi- 
culty. Indeed, these high densities are required in order to 
achieve good convergence of the formulation. As yet, no 
satisfactory answer for this behavior has been found. 

Gi -youty other parameter available to affect the conver- 
sence of the solution is By the wire-to-boundary distance. 
As mentioned in section IITI.8.2., this parameter can be ex- 
pected to be critical in achieving a good field approximation. 
Computed results show this to be true. 

A great deal of effort was expended in looking for errors 
in the formulation before the solution behavior with respect 
to BL and Ny Tow wMiveamermecltated. [Lt 1S still possible, no 
matter how unlikely, that some error still exists. If an 
error is present in the formulation, or code, it's effect is 
Guat]? peculiar. By comparison with other FEM problems, the 


possible errcr causes two to three times the "normal" node 
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density to be needed, but otherwise seems to have no effect. 
Memes, the results are satisfactory in all respects save 
node density required for convergence. 

The primary comparison case chosen was a one-half wave- 
length wire with a radius of 1/50 wavelength. This case was 
a compromise forced by excessive required computing for longer 
Wires and questionable integral equation comparison results 
With shorter wires. These other cases were looked at, and 
some example results are included here. Their investigation, 
however, has not been exhaustively completed due to lack of 


teenie: 


A. INTEGRAL EQUATION COMPARISON 

The basis used for judging the convergence of the FFM 
computation was a standard and successful type of integral 
equation formulation, Hallen's equation for scattering by 
a thin conductive cylinder. This method assumes a Z directed 
current density, and finds this wire current by a simultaneous 
Solution of two expressions for the vector potential, A. One 


equation is just equation (33). The second is a different 


form of equation (20) 
2 2- _ =S;— 
eet KoA = jwe E (r) (47) 


where E-(r) is the scattered electric field. Simple manipu- 


ieelons yield the integral equation for I(z'), [7] 
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L Z, 
ya ! ! : ! Zi ene) . ' ' 
mG“ ECZ 7’ )VdzZ! = 3 QZ!) eee SS ELD ai 

Q 0 


O 


iC oeme woos KZ = UGZ) (48) 


where A is an impedance per unit length, C and D are con- 
Sees Of integration, and U contains the incident field. 
Equation (48) is coded in subroutine INTEQ. 

Results for the 1/2 x 1/50 wavelength case are shown on 
pages A-1 through A-38 of Appendix A. Ee calculations are 


shown as a solid line, while integral equation computations 


are displayed as the "X" curve. The parameters of each case 
are given below the graph. The magnitude plots are scaled 
in milliamperes along x, and wavelengths along y. The wire 


is to be seen lying along the y axis, and is conveniently 
@emearecated by the zero end values of the integral equation 
curve. The current phase plots have the same y axis units, 
and are scaled in degrees along x. 

Graphs A-l through A-24 show FS convergence with an 
increasing number of nodes. Values of Ny are 20,40,70 and 
100, with three plane wave incidence angles (ALPHA) used in 
each case. In these plots, the plane wave is incident from 
the upper right for ALPHA=30 and 60 degrees. ALPHA=90° is 


broadside incidence. Convergence of the current magnitude 


4UQ 





is plainly seen at all incidence angles. For a=90°, the 
error is 15.5% at the wire center for N_j=20, bucmomily 2.5% for 
NA70. Similarly, the phase error decreases from ice 
Comparing the Va and N= 100 cases, it is seen that the 

extra nodes make very little difference in the quality of the 
solution. Convergence has been essentially accomplished at 
Nj=70; and the additional memory and computing time for N 7100 
are not really necessary. The convergence is less satisfying 

for incidence angles other than 90°. At ALPHA=30-, for example, 
the magnitude is not skewed quite enough, while the phase error 
increases toward the wire bottom. This situation does not 
improve as Ae increases. 

Graphs A-25 through A-38 show FFM convergence as a function 
one Bo ae ae calculations are made with N 728. The behavior 
shown in these plots is characteristic of all cases of wire 
kemeth and radius. For BF much too large, the magnitude is 
too low and the phase is off by as much as g0°. As By decreases, 
the magnitude rises, overshoots, and then settles back towards 
the correct, convergent value. The phase moves smoothly down 
to the correct value. Notice that the Be curves are nearly 
feenmiveameror DR=DZ/32 and DR=DZ/64. This implies that 
convergence has been reached by the time DR has decreased to 
Daasc, or B =.0005X.. For this case the convergence point of 


L 


BF 1s actually about BF 


Results for several other cases of wire length and radius 


=.001A_, as shown in A-1 through A-25. 


are shown in Appendix B. The one and two wavelength wire 


ud 


! 








calculations show that more nodes are required, assuming that 


B.=.001X. is still adequate since A, has not changed. If 70 


L L 
nodes are required for a good answer at Z aA fe: then for 

two wavelengths as many as 300 nodes might be necessary. 
Similarly, well over 100 nodes could be needed for the one 
wavelength wire. The one wavelength case indicates that 
this is probably true. The ee SOUeUne tor net Nj zL00 1s percep- 
tibly better than the Nee 0 solution, especially for off-axis 
incidence angles. 

Appendix C contains graphs for two different cases: 


A, =1/100i. while Z =A of 2 or Zp =A G/ 4. Graphs C-l1 through C-4 


L 


again show convergence value is no longer B =.001\., but is 


L 
smaller than B =.0005A.. 


L 

Plots C-7 through C-10 show the quarter-wavelength wire 
case. It appears from the FFM curves that convergence has 
been almost reached for Nate Tieseeatkou into question 
the accuracy of the integral equation solution. This possible 
inaccuracy should not be too surprising, however, since the 
wire is relatively thick for it's length, and the integral 
equation does much better for very thin wires. 

An "error isolation" check was made on the ee system by 
testing the feedback path. A fictitious sinusoidal wire 
current was generated and its radiated field calculated via 
ho=T hi This field was compared to that given in reference 
Mew the Comparison is shown in plots C=-11 and C-12. The 
two field calculations are equal, to the fourth significant 
digit. This indicates that the feedback path formulation and 


program coding is correct. 
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B. ERROR ANALYSIS 
Consideration of the near field behavior reveals the reason 


for the critical dependence of the solution on B The magne- 


1 
tic field magnitude in the very near field is dominated by a 
WR term, as can be deduced from equation (37). This is 
depicted in Figure 10. Also shown are two dashed lines which 
are the linear basis function expansion of the field. It is 
Clear from the figure that the approximation of the field 


deteriorates as AD decreases, if 5 Fo 1s held constant. Thus 


it is seen why B =. 001A. is adequate for A, =.02A. while 


L 


B, =. 0005. to obtain the best convergence with AL 


1s noted, however, that even at this best convergence value 


= Oy) eee 
O 


oF By the fr solution 1s not as close to the integral equation 
solution as in the A.=.02i_, B 
L O L 


also be explained by consideration of the basis function 


=. 001). case. This fact can 


approximation's behavior. While decreasing B, improves the 


L 
field approximation very near the wire, the approximation 
becomes worse in the neighborhood of R=B,. 
It thus appears that there is a limit to the usefulness 
Semenis LOormulation of the fe solution to the wire-scattering 
Pegpiem. The primary difficulty lies in the inaccuracy of 
the linear basis function expansion used in the FEM calculation. 
A much better basis function may contain a term proportional 


to wR This would greatly reduce the sensitivity of the 


Semcon to BL: 
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Eaoure 10 


There is another source of error which is inherent in 
eeeserc Of Computer solution: roundoff. Simple calculations 
Will be affected in the sixth or seventh significant digit by 
the use of single-precision arithmetic on the IBM-360. If 
only a few operations were required, this error could be 
ignored. But when solving large systems of equations many 
thousands of operations are required, and accuracy can be 
severely degraded. 

A rough guide to the degree of trust to be placed ina 
solution can be had in the condition number of an inverted 
matrix. This number provides an upper bound on error propa- 


gation in a system due to the inaccuracy of finite-precision 


rar 








Meeenmcocle. FOr example, in the system X = nee a condition 
momeer of 1000 for en would cause errors in the third signi- 
egies it Omer (© appear as errors in the third, second or 
eee irst Significant digit of xX. 

This condition number is a measure of the similarity of 
the system equations. In the wire-scatterer case, as the 
solution points are placed closer together, the field values 
at these nodes become more indistinguishable. Hence the 
equations relating the field values at neighboring nodes 
become more nearly the same. The system solution, then, 
becomes sensitive to errors in the less significant digits. 

Condition numbers of system matrices for the cases dis- 
cussed here panies from aeeur 250 to 1000. Cases with fewer 


nodes and larger values of B, have lower condition numbers, 


L 


while as Ne. increases and/or B, decreases, the condition 


L 
number increases rapidly. This is another reason why a good 
solution is more difficult to achieve for very thin or long 
wires. The use of a basis function more suited to the field 


behavior would allow larger values of B, to be employed. This 


L 
would lower the system condition number significantly, and 
improve the solution's accuracy and trustworthiness. Other 


than this and keeping nee 


as small as possible, there is 
nothing that can be done to avoid the problems reflected in 
the condition number. 

Even if the errors mentioned above could be eliminated, 


the results of both the F? and integral equation calculations 
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would still not represent the fields with complete accuracy. 
Baewassumptions made in both formulations are not strictly 
true for any wire of finite thickness. Current flows on 
the surface of the wire, not in the center, as assumed in 
the vector potential formulation. 

Pield phase variation from side to side of the wire is not 
zero, so there is anaZimuthal component of current as well as the 
Meeeected current assumed. And because this aaa 
current exists, the current is not zero at the wire ends. 

But, since essentially the same types of approximations are 
made for both solution methods, it is not possible to ascer- 


tain, from the comparisons made here, the degree to which 


either solution represents the actual fields. 
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V. CONCLUSIONS 


As a result of this investigation, the field-feedback 
formulation has been shown to be a sound technique for the 
computation of electromagnetic scattering problems. It has 
been demonstrated that PF calculations compare very favorably 
with results obtained by the conventional integral equation 
Pommmleation. The success of this initial development is not 
unqualified, however. The two deficiencies previously dis- 
cussed, basis function unsuitability and excessive number of 
required nodes, could present formidable difficulties in 
problems where a priori knowledge of field behavior is not 
available, or where bodies large compared to a wavelength are 
eed. 

The discrepancy between required and expected node 
density 1s particularly troublesome, as this fact is as yet 
unexplained. If high node densities are actually required 
by the method and are not just the result of some formulational 
error, then the attendant large matrices will limit the use- 
fulness of eS computations because of the great amount of com- 
puter memory required. 

In any case, the flexibility of the finite-element method 
when used in the coupled-azimuthal potential formulation 
should allow numerous heretofore difficult problems to be 


attacked. The ease with which mixed boundary conditions can 


7 





be met for complicated inhomogeneous scattering bodies, and 
the accuracy and simplicity of the vector potential formu- 
lation promise a wide applicability for the field feedback 


formulation. 
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